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Abstract 

A semi-analytical dynamical mean-field approximation (DMA) has been de- 
veloped for large but finite iV-unit active rotator (AR) networks subject to 
individual white noises. Assuming weak noises and the Gaussian distribu- 
tion of state variables, we have derived equations of motions for moments 
of local and global variables up to the infinite order. In DMA, the original 
A-dimensional stochastic differential equations (DEs) are replaced by three- 
dimensional deterministic DEs while the conventional moment method yields 
(1/2)N(N + 3) deterministic DEs for moments of local variables. We have 
discussed the characters of the stationary state, the time-periodic state and 
the random, disordered state, which are realized in excitable AR networks 
depending on the model parameters. It has been demonstrated that although 
fluctuations of global variable vary as l/y/N when N is increased, those of 
local variables remain finite even for N — > oo. Results calculated with the use 
of our DMA are compared to those obtained by direct simulations and by the 
Fokker-Planck equation which is applicable to the N = oo AR model. The 
advantage and disadvantage of DMA are also discussed. 
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I. INTRODUCTION 



It has been shown that noises play intrigue and essential roles in non-linear systems. 
Some examples are the noise- induced transition and the stochastic resonance 0. Effects 
of noises on the dynamics of a variety of real systems such as the Josephson junction ||, 
chemical reaction ||, charge density of states || and neural networks ||, have been studied 
with the use of the coupled phase model and its variants. These model are described by 
stochastic, nonlinear differential equations (DEs), which have been solved by simulations 
or analytical methods such as the Fokker-Planck equation (FPE) and the moment method. 
In order to make our discussion concrete, we hereafter pay our attention to the active 
rotator (AR) model, which was first studied by Shinomoto and Kuramoto ||. They 
obtained the FPE for the infinite-dimensional (N = oo) AR model, discussing the phase 
transition between the stationary state and the time-periodic state. Responses of excitable 
AR models to subthreshold inputs have been investigated ||- [|12|| . In recent years effects of 
additive and/or multiplicative noises in AR models have been studied by using FPE ||13|| - 
||16|| . Although FPE method is a powerful tool for solving stochastic DEs, its application is 
limited to a single or infinite system subject to white noises. Kurrer and Schulten proposed a 
moment method, expanding FPE for the N = oo model in a Taylor series around the center 
of distribution [T^ . Rodriguez and Tuckwell (RT) adopted a different moment method in 
which the original stochastic DEs are expanded in a series of fluctuations around means 
of state variables [TS|]- [21]. RT's original moment method was first applied to Fitzhugh- 



Nagumo (FN) neuron model |TS| |22] and then FN neuron networks |[L8| . RT's moment 



method takes into account means, variances and covaricance of local variables, then the 
number of DEs for A- unit FN neuron networks is N eq = N(2N + 3), which is, for example, 
20300 for iV = 100. When we apply RT's moment method to the iV-unit AR network under 
consideration, we get N eq = (1/2)N(N + 3), which is 5150 for iV = 100. This exponential 
increase in the number of DEs prevents us from calculations for a system with a realistic 
size. 

Quite recently the present author has proposed an alternative moment method, which 
is hereafter referred to as a dynamical mean- field approximation (DMA) P3|. In DMA we 
take into account means, variances and covariances of local and global variables, replacing 
the original 2 N- dimensional DEs of the iV-unit FN model by eight-dimensional DEs in- 
dependently of N. We have investigated the N dependence of the spike timing precision, 
which has been shown to be improved by increasing N. The purpose of the present paper 
is to apply DMA to coupled AR networks to discuss their dynamics. Although in our pre- 



vious paper | 23fl , we included up to the forth-order moments, we will, in this paper, take 
into account up to infinite-order moments, which are expressed as a product of second-order 
moments with the use of the Gaussian assumption for the distribution of state variables. 
The number of equations for iV-unit AR networks becomes N eq = 3 in DMA, which is much 
smaller than N eq = (1/2)N(N + 3) in the conventional moment method. 

The paper is organized as follows: In Sec. II, we develop a DMA theory for an ensemble 
of N systems, obtaining equations of motions of moments of local and global variables. 
In Sec. Ill, the phase diagram showing various states in coupled AR models is discussed. 
Discussions are given in Sec. IV, where we compare our DMA with the conventional moment 
method, showing that the former may be alternatively derived from the latter with a proper 
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reduction of the number of variables with the mean-field approximation. The final Sec. V 
is devoted to conclusions. 



II. DYNAMICAL MEAN-FIELD APPROXIMATION 

A. Basic formulation 

We assume an ensemble consisting of coupled TV-unit phase models subject to white 
noises, whose dynamics of the phase 0« (mod 2tt) of the ith system is described by nonlinear 
DEs given by 

d -^ = m) + ^T.G^-<Pi)+ut), (< = i-iv) a) 

where explicit forms of F(x) and G(x) will be given shortly [Eq. (14)], w denotes the 
coupling, and £j(t) is the independent Gaussian white noise with < ^(t) >= and < 
£i(t) £j(t') >= 2D 5ij S(t — t') : the bracket < • > denoting the average over stochastic 
random variables [see Eq. (4)]. 

We will express these nonlinear DEs by moments of local and global variables of the 
ensemble. The global variable is defined by 

*(*) = ^ E Mt), (2) 



and their averages by 



where the bracket <> denotes 



fi(t) = < m >, (3) 



<G(0)> =/•••/ d<t>G(4>,t)p(<l>,t), (4) 

p(4>,t) being the probability distribution function (pdf) for the iV-dimensional stochastic 
variable </>=(0i, 4>n)- 

We express DEs given by Eq. (1) in terms of the deviations from their average given by 

sut) = Ht)-Kt), (5) 

to get variances between local and global variables, given by (the argument t is hereafter 
neglected) 

7 = ^ E < ^ >> (6) 

i 

p = <^ 2 >=^E E<<^<%>> (7) 

i 3 

where <5<£> = <£>(£) — /i(t). It is noted that 7 expresses the spatial average of fluctuations in 
local variables of 0, while p denotes fluctuations in a global variable of $. 
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We assume that the noise intensity D is weak. This allows us to expand the right-hand 
side of Eq. (11) around the average p, to get 



rifh 00 i;i 00 City 

^ = E V W + 1 S E - ^ + ^ (8) 



from which we get 

dp 1 x-^ c?0i 
— = — > < — — > 
dt N . dt 

% 

^EE^r <(^) £ >+^EEEV < ^-^) >> ® 

where F® = F^{p) and G® = G®(0). 

We furthermore assume that pdf of state variables takes Gaussian form. Numerical 
simulations have shown that for weak noises, the distribution of <p(t) in a single AR system 
nearly obeys the Gaussian distribution, although for strong noises, the distribution of 4>(t) 
deviates from the Gaussian |fLl |. Similar behavior of the Gaussian distribution of state 
variables has been reported in FN and Ho dgkin- Huxley neuron models p4 |. When 
we adopt the Gaussian decoupling approximation, averages higher than the second-order 
moments in Eq. (9) may be expressed in terms of the second moments given by 

< 5(pi > = n *™ < ^4>k^<Pm >, for even £, 

all parings 

= 0, for odd £, (10) 

where the summation is performed for all (£ — 1)(£ — 3). ...3 ■ 1 combinations. After some 
manipulations with the use of the approximations mentioned above, we get equations of 
motions for p, 7 and p given by (for details see Appendix A) 

ai n=0 U - Z n=0 U - 

4 = 2 7 E — ^(|) n + - p) E —^(7 - PY + 2D, (12) 
at n=Q 11. z n=Q 11. 

dp ^F( 2n+1 ) 7vn 2D . . 

The coupled AR network is given by 

^^ = c -a«n(x) + ^X;«n(^-0 i )+6(t), (< = 1 - JV) (14) 

with F(x) = c — a sin(x) and G(x) = sin(x) in Eq. (1). In Eq. (14) c (> 0) stands for the 
intrinsic frequency and a the intensity of the pinning force introduced such that for c < a, 
the system mimics the stochastic limit cycle or excitable elements. The model with c = 
stands for the equilibrium planar model. The case of a = corresponds to a usual phase 
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model |H p5| . For w = and ^ = 0, the AR system locates at the stationary point given by 
(pi = 0* = arcsin(c/a). When noises are introduced, the system shows the intrigue behavior. 
Substituting F(x) and G(x) to Eqs. (11)-(13), we obtain DEs for p, 7 and p given by 

= c _ a s in{p) exp(~), (15) 
d'y 7 

— = —2a 7 cos(p) exp(— — ) — 2w( r y — p) exp[— (7 — p)] + 2D, (16) 

^ = -2a p cos(p) exp(~) + (17) 

The original iV-dimensional stochastic DEs given by Eq. (1) are transformed to three- 
dimensional deterministic DEs, which show much variety depending on model parameters 
such as a, c, w, D and N. 

We note that the noise contribution is 2D in Eq. (12) while that is 2D/N in Eq. (13). 
It is easy to see that 

p = j/N, (for w/D -> 0) (18) 

= 7. (for D/w -> 0) (19) 

Equation (18) agrees with the central-limit theorem. In the limit of N = 00, we get p = 0. 
On the contrary, in the limit of N = 1, we have p = 7. 



B. Various quantities 



Distribution of local variables 

Adopting the mean-field approximation, we get < (pi >~ (l/N)J2k < 0fc >= A* an d 
< >~ (1/-/V) < 50| >= 7. Then the distribution for the variable (pi is given by 

P(«^(i)^), (20) 

where 0(x) is the normal distribution function given by 

m = 4=exp(-^). (21) 

V 27T Z 

The probability given by Eq. (20) depends on the time because of the time dependence of 
pit) and 7(t). 

Distribution of global variables 

Mean and variance of global variables $ are given by < $ >= p and < 5$ 2 >= p, 
respectively. We get the distribution for the global variable $ given by 

PW ^ ( ^ m ^ y (22) 
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Averaged frequency 

The averaged frequency v is denned by 

" = l 7V(]v^T)?? T -l"- < 23) 

with 

T ik = Uk+i — Uk, (24) 

tik = {t I fa(t) = 9;fa>0;t> t lk _ x + r r }, (25) 

where Nfi stands for the number of firings of a given rotator i, (pi = dfa/dt, T oik the interspike 
interval (ISI) of output signals, t ik the fcth firing time, and 9 (= 2 7r) and r r (= 5) are the 
threshold level and the refractory period, respectively. When there is no firings, we set 
v = 0. 

Order parameters 

The order parameter ( and its fluctuation 5( are defined by 



C=I*(*H (26) 

5( = VUWT - C 2 , (27) 

with 

z (t) = ex v[i fa(t)\ =< ea;p[i >, (28) 

where the overline denotes the temporal avarage. By expanding z(t) in a series of 8 fa around 
\i and adopting the Gaussian decoupling approximation given by Eq. (10), we get 

z(t) = exp[i v(t)] £ = exp[< /!(*) - ^] (29) 

n n. z z 



Synchronization ratio 

The synchronization ratio a is defined by 



a = s(t), (30) 

with 

For the completely synchronous (asynchronous) state, both £ and cr become 1 (0). It is 
noted, however, that while ( depends on 7, a is a function of (p/7 — l/N) = iV _1 < 
5 fa 8(f) j > /J2i < fitfii >! the ratio of the inter-AR correlation to the intra- AR correlation. 
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Before discussing calculated results with the use of DMA, it is worth to mention the 
calculation of Kurrer and Schulten [17|. They intended to solve the FPE for the N = oo 
AR model given by |7| 

- n{<f>, t) = -q7[F{4>, t) - w J d(j>' sin{<t> - (j)')n{<t>', t)] n(<j>, t) + D— n{<f>, t), (32) 
where the density probability n((j), t) is defined by 

n(&*) = i£<W -&(*)), ( 33 ) 

i 

with the periodic condition: n(<f> + 2ir,t) = n(<j),t) and the normalization condition: 
/ dcf) n(4>,t) = 1. Kurrer and Schulten |T7| expanded F(<ft,t) in a Taylor series around 
the center of distribution, assuming the Gaussian form for n((f>, t) given by 

n{M = ^^ e JJi^f\, (34) 

where the mean u(t) and variance v(t) obey DEs given by 

— = c — a sin(u) exp(— — ), (35) 

— = —2a v cos(u) exp(——) — 2w v exp(-v) + 2D. (36) 

Equations (35) and (36) resemble our Eqs. (15)-(17) if we read u — > fi and v — > 7. Actually, 
Eqs. (35) and (36) are equvalent to Eqs. (15) and (16) in the case of p = 0, which is realized 
in the limit of N = 00. 



III. CALCULATED RESULTS 

DMA equations given by Eqs. (15)-(17) have been solved by the forth-order Runge- 
Kutta method with a time step of 0.01, the initial conditions being /x(0) = 7(0) = p(0) = 0. 
Calculations have been performed for < t < 1000 (100 000 steps) and results for t < 100 
are discarded. Simulations of directly solving Eq. (1) have been made by the forth-order 
Runge-Kutta method with a time step of 0.01, the initial conditions being 0^(0) = (i — 1 
to N). The number of trials for a given set of parameters in our simulations is hundred 
otherwise noticed. We have solved also FPE given by Eq. (32), which is valid for the 
N = 00 AR model. We first Fourier transform FPE with the first 30 modes after Ref. [|7J. 
A set of 61 ordinary DEs has been solved by the Runge-Kutta method. 



A. Phase diagram for various types of solutions 

By solving Eqs. (15)-(17), we get the stationary state and the non-stationary state: in 
the former state the variables are time independent while in the latter state they are time 
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dependent. The equilibrium values of p, 7 and p in the stationary state are given by (we 
set c = 1 hereafter) 

p = arcsin[{—) exp(^)], (37) 
7 = D + wpea:p[-(7-p)] g 

y/a 2 exp(— 7) — 1 + u> exp[— (7 — p)] 

, = , J/iV . (39) 
'a 2 exp(— 7) — 1 



The stationary state where Eqs. (37)-(39) are satisfied, is hereafter referred to as the S state. 
Kurrer and Schulten |T7] pointed out that the non-stationary state may be classified into 
the time periodic (P) state and the random, disordered (R) state. DMA also yields three 
types of the S, P and R states characterized by the quantities of (, 5(, v and a introduced in 
Sec. IIB, result being summarized in Table 1. In the S state, 5( and v are vanishing while 
( (~ 1) and a are finite. In the P state, all quantities are finite. In contrast, in the R state, 
all quantities except v vanish. 

Boundaries between these three states depend on a, w, D and N. Figure 1 expresses the 
D — a phase diagram showing the boundaries between these states in coupled AR models 
calculated with the use of DMA for iV = 10, 100 and N = 00. The gradient of the boundary 
between the stationary (S) state and non-stationary (P+R) states is decreased as increasing 
the value of w and/or of N. The difference between boundaries for N = 10 and N = 00 
with w — 0.1 is very small: the effect of N becomes more significant for a stronger coupling. 
The critical a value, a c , above which the S state exists, is given by 

a c - 1 ~ [ Cl - c 2 w (1 - 1)] D, (40) 

where C\ = 2.25 and c 2 = 1.75. In contrast, the critical value of ad for the boundary between 
the P and R states for w — 1.0 is given by 

a d -l^-(d 1 + ^)[D-(d 3 + ^)], (41) 

where d± = 5.36, d 2 = 257, d 3 = 0.265 and d± = 0.9. 

The behavior of the solutions of DMA in the S, P and R states when D and/or a values 
are changed, is shown in Figs. 2-5. We will first mention the calculations of DMA in 
the three states. Figure 2(a) and 2(b) express the D dependence of (, 5(, v and a for 
a representative set of parameters of a = 1.05, w — 1.0 and iV = 100, showing that the 
network is in the S state for D < 0.082, in the P state for 0.082 < D < 0.273, and in the R 
state for D > 0.273. In contrast, Fig. 3(a) and 3(b) express (, 5(, v and a as a function of 
a for a set of parameters oi D — 0.1, w — 1.0 and iV = 100, for which the networks is in the 
S state for a > 1.06 and in the P state for a < 1.06. 

Equations (37)-(39) for small D and w in the S state yield 

p = arcsin(-) + d\ D + .., (42) 
a 



S 



D 



[l + d 2 D-d 3 (l 



1 



) w + d 4 w 2 + ..], 



(43) 



7 = 



(D/N) 



N 



P 



[l + d 2 D + ..,\, 



(44) 



leading to 




(45) 



1 d s w 



(46) 



N (l + d 2 D) 



where d 1 = l/2(a 2 - 1), d 2 = a 2 /2(a 2 - 1) 3/2 , d 3 = \j\Ja 2 - 1, and d 4 (> 0) is a complex 
function of D, w and N. 

Solid curves in Figs. 4(a) and 4(b) show distributions of local [P((f>i(t))] and global 
variables [P($(t))], respectively, in DMA for a = 1.05, w = 1.0, D = 0.05. They are 
obtained by Eqs. (20) and (22) with p = 1.339, 7 = 0.04354 and p = 0.00212. 

Figures 2(a) and 2(b) show that when the noise intensity is increased and crosses the 
value of 0.082, the AR network begins correlated firings. This implies the appearance of 
the P state, where 5(, v and a are continuously changed. Solid curves in Fig. 4(c) ad 4(d) 
express P(<f>i(t)) and P(<£>(t)) for D = 0.10, respectively, which are given by Eqs. (20) and 
(22) with p = 1.497, 7 = 0.11022 and p = 0.009443. The time evolution of the probability 
of P(<f)(t)) calculated in DMA for D = 0.10 in the P state is shown in Fig. 5(a), which is 
oscillating with the period of about 40. It is noted that not only the position of P(<f)(t)) but 
also its width change as a function of t. For example, we get p = 6.151 and 7 = 1.711 at 
t = 120 while p = 1.497 and 7 = 0.11022 at t = 100. 

When the D value is rmore increased, the AR network fires abundantly and irregularly, 
which suggests the appearance of the R state. The solution of Eqs. (37)-(39) in the R state 
for a large t is given by 



which lead to vanishing (, 5( and a except v. Figure 2(a) and 2(b) show that ( and S( 
suddenly vanish at D = 0.273 with no hysteresis. Solid curves in Fig. 4(e) ad 4(f) express 
P((f>i(t)) and P($(t)) for D = 0.30, respectively, which are given by Eqs. (20) and (22) with 
p = 1.339, 7 = 19.9 and p = 0.199. 

In the following, results of DMA will be compared with those of simulations and FPE. 
Dashed curves in Figs. 2(a) and 2(b) show the results of simulations for the D dependence 
of (, 5(, v and o. The agreement of ( in DMA with that in simulations is good for S 
and P states. However, it is not good in the R state, where ( vanishes in DMA but not 
in simulations. This is expected to be due to deviations of the state-variable distributions 
from the Gaussian form. When D is more increased in the R state, our simulations yields 
a gradual decrease in C, which is 0.454, 0.276, 0.188, 0.139 and 0.110 for D = 1.0, 2.0, 3.0, 
4.0 and 5.0, respectively, with a = 1.05, w = 1.0 and iV = 100. We note in Fig. 2(a) 



p ~ ct, 
7 ~ 2Dt, 



(47) 
(48) 

(49) 
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that 5( of simulations is about ten times smaller than that of DMA. This is clearly seen in 
Fig. 6(a) where we plot the time evolution of | z(t) | obtained by DMA and simulations for 
a = 1.05, w = 1.0, D = 0.10 and N = 100. The former has larger temporal fluctuations than 
the latter although both yield similar averaged values of ( = \ z(t) |. In contrast, Fig. 6(b) 
shows the time dependence of s(t) calculated by DMA and simulations for a = 1.05, w = 1.0, 
D = 0.10 and N = 100. Again our DMA yields larger fluctuations in s(t) than simulations 
although both methods lead to similar averaged values of a = s(t). When comparing K(a) 
with K(b), we notice that the time dependence of | z(t) | is not the same as that of s(t). 
This is because z(t) is a function of 7 [Eq. (29)] while s(t) is a function of the ratio of p/7 
[Eq. (31)]. Dashed curves in Figs. 3(a) and 3(b) show the a dependence of £, 5(, v and 
cr obtained by simulations. P{4>) and -P($) obtained by simulations are plotted by dashed 
curves in Fig. 4(a)-4(f). Dotted curves in Fig. 4(a), 4(c) and 4(e) denote n(0) obtained 
by FPE for the N = 00 AR model. Figure 5(b) express the time evolution of n((j),t) in the 
P state obtained by FPE. From a comparison between the results of DMA and simulations 
(and FPE) mentioned above, we note that DMA is good for the S state, in fairly good for 
the P state in the qualitative sense, but not good for the R state. 

For a comparison, we show by the dotted curve in Fig. 1, the boundary obtained by 
Shimokawa and Kuramoto (SK) with the use of the FPE for w = 1.0 and N = 00 M. The 
ordered P state where 5( 7^ and v 7^ is reported to exist in the triangle region enclosed by 
the dotted curve and the horizontal axis. The P state obtained by SK is nearly in agreement 
with our P state. In SK, states besides the P state are regarded as the stationary state where 
dn((p,t)/dt = 00. On the contrary, Kurrer and Schulten (KS) distinguished the R state 
from the P state, both of which are non-stationary [y 7^ 0) [|H]]. The results of SK and KS 
are for the N = 00 AR model. Figures 7(a) and 7(b) show the D dependence of (, 5( and 
v for a set of parameters of a = 1.05, w = 1.0 and N = 10 4 , which are the same as in Figs. 
5(a) and 5(b) except N. In order to simulate the N = 00 limit, the N value in Figs. 7(a) 
and 7(b) is chosen to be very large but finite because s(t) given by Eq. (31) is not properly 
defined in this limit. Results for N = 00 in Figs. 1 and 7 should be compared with those 
obtained by SK and KS. As was pointed in Sec. IIB, DEs of DMA given by Eqs. (15)-(17) 
in the limit of N = 00 agree with those of KS given by Eqs. (35)- (36). Nethertheless, our 
D — a phase diagram for iV = 00 in Fig. 1 does not agree with that of KS. For example, KS 
obtained the critical values given by (17[ 



a c — 1 = , (between S and P + R) (50) 

2 w 

— = 0.736, (between P and R) (51) 

w 

which do not agree with our expressions given by Eqs. (40) and (41) for N = 00. 



B. Cluster-size (N) dependence in the S state 

Since one of the advantages of DMA is that we can discuss the finite-iV property of 
coupled AR networks, we have made more detailed calculations of the iV dependence of the 
quantities in the S state. Figures 8(a) and 8(b) show the log- log plot of the N dependence of 
7, p and a, results for N = 10 and 20 being for 500 trials. Solid curves in Fig. 8(a) express 
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the result of DMA for a set of parameters of a = 1.05, w = 1.0 and D = 0.05, whereas circles, 
squares and triangles denote those of simulations. For this set of parameters, the P state is 
realized for N < 9. We note that as iV is decreased from above and approaches to the S-P 
boundary, fluctuations of 7 and p are increased (and a is also increased). Similar behavior 
is observed for a different set of parameters. Figure 8(b) shows the results for a = 1.20, 
w = 1.0 and D = 0.10. With this set of parameters, we get the S state for N > 2 and the P 
state for N — 1. Figures 8(a) and 8(b) show that as increasing N, p is much decreased but 
7 shows only a weak N dependence. We note that a is decreased as increaing N whereas 
( = exp(— 7/2) shows little N dependence. We should stress that although fluctuations of 
global variables is inversely decreased as p oc 1/N consistent with the central-limit theorem, 
those of local variables remain finite even for N — > 00. 



IV. DISCUSSIONS 

We have proposed DMA theory for stochastic, nonlinear networks like coupled AR mod- 
els, taking into account means, variances and covariances of local and global variables. It is 
worth to compare DMA with the conventional moment method in which means, variances 
and covariances of local variables are given by 

m { = < fa >, (52) 
dj = < A<& A<f)j >, (53) 

with A(pi = (pi — mi. By using the Gaussian decoupling approximation [Eq. (10)], we get 
(for details see Appendix B) 

dm 00 F (2n) w 00 <7 (2n) 

— = T — C™ + -YY \C kk + C u - 2C ik ] n , (54) 

dC ■ 00 F (2n+1) w 00 (7( 2n+1 ) 

dt ~ h 2" n\ + 5 + N \ ^ 2" n\ 

x [{Ckk + Cu — 2C ik ) n {Cjk — Cij) + {Ckk + Cjj — 2Cjk) n (C ik — CV,-)] + 2D i 5 i j ) (55) 

where F { / ] = F^( mi ) and = G®(0). 

For the AR model, Eqs. (54) and (55) become 

= c - a sin(mi) exp(~Cu), (56) 
dC 00 ( — )) n 

-JL = -a cosfa) £ V^l ^ + ( n )C - + 2D6 ^ 
w 00 ( — l) n 

k n=0 

+(Cjj + Cfcfc — 2Cjk) n {Cik — Cij)}, (57) 
For variances (i = j), Eq. (57) becomes 
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dC 1 

= -2a cos{rrii) C u exp(--Cu) + 2D 

2w 00 (— l) n 

k n=0 

Taking into the symmetry relations: CV,- = Cji, we get the number of DEs in the moment 
method to be N eq = N(N + 3)/2, which is 65, 5150 and 501 500 for N =10, 100 and 1000, 
respectively, while N eq = 3 in our DMA. 

It will be shown that we can derive DMA from the moment method by reducing the 
number of DEs, adopting the mean-field approximation: 

rrii ~ p, (59) 
d ~ 7 n_1 Cjj, (60) 
(C kk + C« - 2C ik ) n ^ 2™- 1 ( 7 - p) n -\C kk + C« - 2C lk ), (i + k) (61) 



with the relations given by 



^=^E m *> ( 62 ) 



N 

7 = ^E C '«» (63) 



p = ^EE^v (64) 

* 3 

DEs for /i, 7 and p are given by Eqs.(ll)-(13) for the general phase model or by Eqs. 
(15)-(17) for the AR model. 

It is possible to regard DMA as the single-site self- consistent theory. Let us assume a 
configuration in which a single nonlinear system % is embedded in an effective medium whose 
effect is realized by a given system % as its effective external input through the coupling w. 
Assuming m 8 = p, we replace quantities in coupling terms of Eqs. (57) and (58) by effective 
quantities of p, 7 and p, to get 

dm, t -F( 2 ") w -G( 2n ) 

*- = iv ?S— <^-'>" (65) 

oY7 00 F( 2n+1 ) 7/; 00 (7( 2n+1 ) 

^ = e w (G & + c «>" c « + ^ ? 5 - < 66 > 

We should note that and Cjj determined by Eqs. (65) and (66) are functions of p, 7 
and p. In order to self-consistently determine them, we impose the self-consistent conditions 
given by 

p = m h (67) 
7 = Cu, (68) 

P=^E^r (69) 

j 
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Note that Eqs. (67)- (69) are assumed to hold independently of i and that rrii and CV,- in 
their right-hand sides are functions of 7 and p. The condition given by Eqs. (65)-(69) 
with the mean-field approximation given by Eq. (59)-(61) yields DEs for 7 and p which are 
again given by Eqs. (11)-(13). The self-consistent condition given by Eq. (67)-(69), which 
assumes that the quantities averaged at a given site are the same as those of the effective 



medium, is common in mean-field theories such as the Weiss theory for magnetism p6[ and 
the coherent-potential approximation for random alloys |27|| . 

By using DMA, we have investigated the response of the excitable, coupled AR networks 
to an applied periodic spike, by adding to the right-hand side of Eq. (15), the input term 
given by 

hnif) =9, f° r m T p < t < m T p + T w (m: integer) 

= 0, otherwise (70) 

where g denotes the magnitude, and T p (=50) and T w (=5) stand for the period and the 
duration of spikes, respectively. We get the critical value of g c = 0.159 below which there 
are no firings for D — 0. Figure 9(a) shows the distribution of ISI, T Q , of output signals 
defined by Eqs. (24) and (25) as a function of D in the absence of input spikes (g = 0) for 
a = 1.05, w = 1.0 and N = 100. Firings begin at D = 0.082, above which the system is in 
the P state as discussed in Sec. IIIA. Around the P-R transition at D = 0.273, ISIs have 
a small distribution. When the input spike is applied, distributions of ISIs are significantly 
changed. Figure 9(b) shows the distribution when the subthreshold input with g — 0.1 
(< g c ) is applied. Firings occur at D > 0.04 with a help of noises. A flat segment at 
0.04 < D < 0.08 corresponds to a periodic solution locked to input spikes while the others 
show the complex behavior. In contrast, Fig. 9(c) shows the distribution of ISIs for the 
suprathreshold input with g = 0.2. At 0.05 < D < 0.08 in the S state, a new branch with 
35 < T D < 43 appears beside the branch with T Q = 50 locked to inputs. The distribution of 
ISIs in the presence of input spikes has much variety than that in the absence of noises, in 
particular in the P state, where the bifurcation is realized as the noise intensity is changed. 
It is possible to discuss the firing-time accuracy of excitable AR models for an external 



input with the use of DMA |[23]| . The fcth firng time of a given rotator i is defined as the 



time when (f>i(t) crosses the threshold 6 from below [Eq. (25)]: 

t ik = {t I <Pi{t) =0]<Pi>0;t> tik-i + r r }. (71) 

By using the discussion presented in Sec. IIB, the probability Wt when <pi{t) at t is above 
the threshold 9 is given by 

Wt(t) = l-il>fiz£) } (72 ) 

where ip{y) is the error function given by integrating the normal distribution function <p{x) 
from —00 to y [Eq. (21)]. The fraction of a given rotator i emitting output at t is given by 

Z,(i) = ^ 8(A) = «^i)|(^)e(A), (73) 
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where Q(x) = 1 for x > and otherwise, and fi = dfi{t)/dt. When we expand //(£) in Eq. 
(73) around t* where nit* ) = 9, we get 

?(^)4(-^)e(/o, (74) 

with 

SU = (75) 

We note that provides the distribution of firing times, showing that most of firings locate 
in the range given by 

t i G [t* a - Stat, t* + 8^]. (76) 
In the limit of vanishing D, Eq. (74) reduces to 

Z e (t) = 6(t-t:). (77) 
Similarly, we define the kth firing time relevant to the global variable $(£) as 

t gk = {t | Ht) =0;6 l >0;t> t gk ^ + r r }. (78) 
The distribution of firing times t g is given by 

z 'W~«^l^ e W- (79) 

with 

5t 03 = (80) 
Equation (79) shows that most of t og locate in the range given by 

tag G [£* - St g, t* a + St og ]. (81) 

From Eqs. (75) and (80), we get 

t, 



• og 



tot V 7 ViV 



(as w/D -> 0) (82) 



This implies that the firing-time accuracy is improved as the ensemble size is increased even 
when there no couplings among ARs. This is consistent with results reported prevously |28| - 

III- 
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V. CONCLUSIONS 



We have developed DMA, which has been shown to be derived in various ways: equations 
of motions for means, variances and covariances of local and global variables (Sec. IIA), a 
reduction in the number of moments in the moment method, and a single-site self-consistent 
approximation to the moment method (Sec. IV). Our DMA theory, which assumes weak 
noises and the Gaussian distribution of state variables, goes beyond the weak coupling 
because no constraints are imposed on the coupling strength. The advantage of DMA is 
that it can be applied to large but finite-A" nonlinear systems subject not only to white 
noises but also to color noises. This is in contrast with FPE, which is applicable to a single 
or infinite system subject to white noises. The limitation of DMA is the weak noise, for which 
calculated results based on DMA are in fairly good agreement with those obtained by direct 
simulations. When the noise intensity becomes stronger, the state-variable distribution 
more deviates from the Gaussian form, and the agreement of results of DMA with those 
of simulations becomes worse. Nevertheless, our DMA is expected to be meaningful for 
qualitative or semi- quantitative discussion on the properties of coupled nonlinear systems. 
It is possible to regard DEs given by Eqs. (15)-(17) as the mean-filed AR model which 
may show interesting behavior for applied input signals and noises. When we consider an 
ensemble of A-unit systems, each of which is described by a M-variable nonlinear DE, the 
number of the deterministic DEs is N eq = M + M(M + 1) = M(M + 2) independently 
of N in DMA while it is N eq = AM + (1/2)NM(NM + 1) = (1/2)NM(NM + 3) in the 
conventional moment. In the case of M = 2 (as in FN model), for example, DMA leads to 
N eq = 8 while the moment method yields N eg = A(2A + 3), which is 2310, 20 300 and 2 
003 000 for A = 10, 100 and 1000, respectively. These figures clearly show the advantage 
and feasibility of our DMA theory. 

To summarized, the property of excitable AR networks has been discussed with the use of 
DMA. The obtained results are summarized as follows. (1) Depending on model parameters 
of a, w, D and A, AR networks display three types of dynamics (Fig. 1): S, P and R states 
are characterized by the quantities of (, 5(, v and a, as summarized in Table 1. (2) The 
S-P transition is of the (continuous) second-order one while P-R and S-R transitions are of 
the (discontinuous) first-order one with no hysteresis. (3) There are no enhancements in 
order-parameter fluctuations of 5( at the transitions. (4) Fluctuations in local variables (7) 
remain finite even for N = 00 whereas those (p) in global variables varies as p oc 1/A, which 
is consistent with the central-limit theorem. 
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APPENDIX A: DERIVATION OF EQS. (11)-(13) 

When we adopt the Gaussian decoupling approximation given by Eqs. (10), Eq. (9) 
becomes 

+ a^EE E (2^1 ^ < ^ - 6 ^ 2 >n ' ( A1 ) 

where B 2n = (2n — l)(2ra — 3). ...3 • 1. In deriving Eq. (Al), we treat (5<f>j — 5<f>j) as a new 
variable with the Gaussian distribution. Adopting the mean-field approximation given by 

< 8$ > n ~ 7"- 1 < ^ 2 >, (A2) 
< (50, - <50,) 2 > n ~ 2"" 1 (7 - p)"- 1 (< 50 2 > + < 50? > -2 < 50, 50, >). (i ^ j) (A3) 



we get 



r/i/ 00 /?(2n) 00 /nf(2n) 

ai n=0 Z n=0 



which yields Eq. (11). 

From Eqs. (8) and (9), we get 

d5(pi d(pi dji 



dt dt dt ' 



(A5) 



= V- F (2n+1) l"Vt>> , V- F (2n)r 

to (2n + l)! t [ (2n)! 2™n! J 

, E v V r (2n + D W - m 2n+l 

Nyto " (2n+l)! " 

+ ^tnt [_ (2^)! 

With the use of Eq. (A6), the calculation of d^/dt is performed as follows. 



d'j 2 ^ 

9 00 p(2n+l) 9 

4?S(£TT)T ( <^ +2 >^?<*^> 

< >, (A7) 

9 00 p(2n+l) 

9 00 Q(2n+1) 

- W2 E E E (2^1)? 5 - +2 < ^» - >™ +1 • (A8) 
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By using the mean-field approximation given by Eqs. (A2) and (A3) and 

< 5^(5^ - <%) >" +1 ^ (7 - p) n (< 5$ > - < >), (i + j) (A9) 

we get 

J oo p(2n+l) oo fi{2n+\) 

ai n =0 Z ' 4 - n=0 IL - 

which leads to Eq. (12). 

The calculation of dp/t is similarly performed by 

dp 1 c.dSc^j d5(f) i 

Tt = w T,T,<^ + ^>. (aid 

which yields Eq. (13). 

For the AR model given by Eq. (14) with F(x) = 1 — a sin(x) and G(x) = sin(x), we 

get 

F {e \p) = c-asin(p), (£ = 0) 

= (-l) n+1 a sin(p), (£ = 2n> 0) 

= (-l) n+1 a cos(/x), (£ = 2n + l) (A12) 

GW(0) = 0, [l = 1n] 

= (-l) n , (£ = 2n+l) (A13) 

which yield Eqs. (15)-(17). 

APPENDIX B: DERIVATION OF EQS. (54) AND (55) 

The moment method takes into account means, variances and covariances defined by 

mi = < ^ >, (Bl) 
C ij = <A(f> i A<l> j >, (i,j = ltoN) (B2) 

where A0j = 0j — m 8 and Ca denotes variances. By adopting the Gaussian decoupling 
approximation given by Eq. (10), we get 



drrii 
dt 


oo 

= E 

n=0 


_p(2«) 

(2n)! 


dC^ 


oo 

= E 

n=0 


p(2n+l) 


dt 


(2n+ 1)! 



(B3) 



fi 2n+2 (< A$ >" + < A# >») < A(j>iA<j)j > 

+>77 E E ToTTW^ 2 ^ (A0fc - A ^ )2 < A ^ (A0fc - A< ^ > 

+ < (A(f> k - A&f > n < A^-(A0 fe - A0 4 ) >] + 2D S ij: (B4) 

where = F^(m,i) and G^ = G^(0). By a proper re- arrangement, Eqs. (B3) and (B4) 
reduce to Eqs. (54)- (55). 
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type of states 


c 




V 


a 


S 


F 








F 


R 








F 





P 


F 


F 


F 


F 



Table 1 Quantities in the S, R and P states of coupled AR networks: F and denote 
finite and vanishing values, respectively. 
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FIGURES 



FIG. 1. The phase diagram of coupled AR networks, showing the stationary (S) state, the 
time-periodic (P) state and the random, diordered (R) state, calculated by DMA with N = 10 
(thin solid curves), N = 100 (solid curves) and N = oo (dashed curves). The dotted curve denotes 
the boundary obtained by SK (Ref.[7]). Calculations by changing a parameter of D (a) along the 
horizontal (vertical) chain curve, are presented in Fig. 2 (Fig. 3). 

FIG. 2. The D dependence of (a) ( and S(, and (b) u and a, for a = 1.05, w = 1.0 and N = 100: 
solid curves denote results calculated with the use of DMA: circles ((), diamonds (5( x 10), squares 
[y x 10) and triangles (cr) express results obtained by simulations, dashed curves being drawn only 
for a guide of the eye. 

FIG. 3. The a dependence of (a) £ and 5^, and (b) v and a, for D = 0.1, w = 1.0 and N = 100: 
solid curves denote results calculated with the use of DMA: circles (£), diamonds (5( x 10 ), squares 
[y x 10) and triangles (cr) express results obtained by simulations, dashed curves being drawn only 
for a guide of the eye. 

FIG. 4. The distribution of local [P(<f>(t))] and global variables [P(<f>(t))] for D = 0.05 [(a) and 
(b)], D = 0.10 [(c) and (d)] and D = 0.30 [(e) and (f)] (in arbitrary units). Dashed curves express 
simulation results. Dotted curves in (a), (c) and (e) denote the results of FPE for N = oo. 

FIG. 5. (a) The time evolution of P{(f>(t)) calculated by DMA for a = 1.05, w = 1.0, D = 0.10 
and ./V = 100, and (b) the time evolution of n((p,t) calculated by FPE for a = 1.05, w = 1.0, 
D = 0.10 and N = oo (in arbitrary units). 

FIG. 6. The time evolution of (a) | z(t) | and (b) s(t) for a = 1.05, w = 1.0, D = 0.10 and 
iV = 100: solid and dashed curve denotes the results of DMA and simulations, respectively. 

FIG. 7. The D dependence of (a) ( and 5^, and (b) v and a, for a = 1.05, w = 1.0 and 
N = 10000 calculated with the use of DMA. 

FIG. 8. The log-log plot of the N dependence of 7, p and a for (a) a = 1.05 and D = 0.05 and 
(b) a = 1.20 and D = 0.10, with w = 1.0 and N = 100. Solid curves denote results of DMA, and 
Circles (7), squares (p) and triangles (a) express those of simulations, dashed curves being only 
for a guide of the eye. Right vertical scales are for a only. 

FIG. 9. Distributions of output ISIs, T a , as a function of D for (a) g = 0, (b) g = 0.1 and 
g = 0.2, with a = 1.05 and iV = 100. Arrows in (a) denote the S-P and P-R transition points. 
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